######### This produce Table 3 from the Cleaned up data file (CleanUpFile_Experiment). #########

library(reshape)
library(ggplot2)
library(car)
library(foreign)
library(xtable)
library(stargazer)
library(psych)


stderr <- function(x) sqrt(var(x)/length(x))


load("Experiment_Output.RData")

data1$ageR <- (data1$age-18)/71
data1$incomeR <- recode(data1$income, "5=1;15=2;25=3;35=4;45=5;62.5=6;87.5=7;125=8;175=9")
data1$incomeR <- (data1$incomeR-1)/8
data1$educR <-(data1$educ-1)/4
data1$pidR <- recode(data1$pid7, "-3=6;-2=5;-1=4;0=3;1=2;2=1;3=0")
data1$pidR <- data1$pidR/6
data1$female <- 1-data1$male

data1$govt_construct <- (data1$govt_grid_2+data1$govt_grid_3+data1$govt_grid_4 + data1$govt_grid_5+data1$govt_grid_6+data1$govt_grid_8+ data1$govt_grid_1+data1$govt_grid_7)/8
data1$govt_constructR <- (data1$govt_construct-1)/4


data1$vict_construct <- (data1$ht_vict_1 + data1$ht_vict_2+ data1$ht_vict_3+ data1$ht_vict_4+data1$ht_vict_5+ data1$ht_vict_6+ data1$ht_vict_7)/6
data1$vict_constructR <- (data1$vict_construct-1)/4

govt_c <- lm(govt_constructR ~ as.factor(treat), data=data1)
govt_alldem <- lm(govt_constructR ~ as.factor(treat)+ageR+female+incomeR+white+pidR+educR, data=data1)


vict_c <- lm(vict_constructR ~ as.factor(treat), data=data1)
vict_call <- lm(vict_constructR~ as.factor(treat)+pid7+income+educ+male + white + age, data=data1)

data1$behavior_construct <- (data1$behavior_1 + data1$behavior_2 + data1$behavior_3 + data1$behavior_4 + data1$behavior_5 + data1$behavior_6 + data1$behavior_7)/7
data1$behavior_constructR <- (data1$behavior_construct-1)/4

behavior_c <- lm(behavior_constructR~ as.factor(treat), data=data1)
behavior_cALL <- lm(behavior_constructR~ as.factor(treat)+ageR+female+incomeR+white+pidR+educR, data=data1)

data1$ht_govt2 <- recode(data1$ht_govt, "1=5;2=4;3=3;4=2;5=1")
data1$ht_govtR2 <- (data1$ht_govt2-1)/4

ht_govtR <- lm(ht_govtR2 ~ as.factor(treat), data=data1)
ht_govtR_all <- lm(ht_govtR2 ~ as.factor(treat)+ageR+female+incomeR+white+pidR+educR, data=data1)


govt_mean <- mean(data1$govt_constructR)
vict_mean <- mean(data1$vict_constructR)
behavior_mean <- mean(data1$behavior_constructR)

stargazer(behavior_c, govt_c, vict_c, type="latex",
          column.labels = c("Respondent Behavior", "Government Policies", "Victim Assistance"),
          covariate.labels = c("Security", "Labor", "Sex","Immigration","Domestic"))

stargazer(ht_govtR, behavior_c, govt_c2, vict_c,
          type="latex",
          column.labels = c("Respondent Behavior", "Government Policies", "Victim Assistance"),
          covariate.labels = c("Security", "Labor", "Sex","Immigration","Domestic"))


stargazer(ht_govtR, ht_govtR_all, behavior_c, behavior_cALL, govt_c, govt_alldem, 
          type="latex",
          column.labels = c("Priority","", "Respondent Behavior","", "Government Policies", ""),
          covariate.labels = c("Security", "Labor", "Sex","Immigration","Domestic",
                               "Political ID", "Income","Education","Male", "White","Age"))



###Summary Table for Indices


newdata <- data1[c(167,88:94)]
newdata<- na.omit(newdata)

vict <- newdata[,2:8]
vict.pca <- prcomp(vict, center=TRUE, scale.=TRUE)
plot(vict.pca, type = "l")
biplot(vict.pca)


newdata <- data1[c(167,77:83)]
newdata<- na.omit(newdata)

behavior <- newdata[,2:8]
behavior.pca <- prcomp(behavior, center=TRUE, scale.=TRUE)
plot(vict.pca, type = "l")
biplot(vict.pca)


##### 
govt <- data.frame(data1$govt_grid_1,data1$govt_grid_2,data1$govt_grid_3,
                   data1$govt_grid_4,data1$govt_grid_5,data1$govt_grid_6,
                   data1$govt_grid_7,data1$govt_grid_8 )
govtC <- data.frame(data1$govt_grid_2,data1$govt_grid_3,
                    data1$govt_grid_4,data1$govt_grid_5,data1$govt_grid_6,
                    data1$govt_grid_8)
a.govt <- alpha(govt)
a.govtC <- alpha(govtC)

govt <- data.frame(data1$govt_grid_2,data1$govt_grid_3,
                   data1$govt_grid_4,data1$govt_grid_5,data1$govt_grid_6,data1$govt_grid_8 )

alpha(govt)


behavior <- data.frame(data1$behavior_1,data1$behavior_2,data1$behavior_3,data1$behavior_4,
                       data1$behavior_5,data1$behavior_6,data1$behavior_7)
a.behavior <- alpha(behavior)

victim <- data.frame(data1$ht_vict_1,data1$ht_vict_2,data1$ht_vict_3,data1$ht_vict_4,
                     data1$ht_vict_5,data1$ht_vict_6,data1$ht_vict_7)
a.victim<- alpha(victim)

report <- data.frame(
  obs= c(length(data1$behavior_constructR), length(data1$govt_constructR), length(data1$vict_constructR)),
  means = c(behavior_mean, govt_mean, vict_mean),
  stdv = c(sd(data1$behavior_constructR), sd(data1$govt_constructR), sd(data1$vict_constructR)),
  mins = c(min(data1$behavior_constructR), min(data1$govt_constructR), min(data1$vict_constructR)),
  maxs = c(max(data1$behavior_constructR), max(data1$govt_constructR), max(data1$vict_constructR)),
  alphas = rbind(a.behavior$total[2],a.govt$total[2],a.victim$total[2])
)
xtable(report)